

% initial guess

c1 = 14/15;
c2 = 0.94;

% tolerance
tol = 0.01;

x = load('tableData.dat');

Hnumerical = x(:,2) ./ 2;

o = optimset('MaxFunEvals',4000,'MaxIter',4000,'Display','off');

s1 = 10;
while(s1 > tol)
    fun1 = @(x, xdata) fitFunctionForHVariable(x,xdata,c1);
    fun2 = @(x, xdata) fitFunctionForHVariable(x,xdata,c2);
    N = [2:65]';
    x11 = lsqcurvefit(fun1,[0,0,0,0,0,0],N,Hnumerical(1:64),[],[],o);
    x12 = lsqcurvefit(fun2,[0,0,0,0,0,0],N,Hnumerical(1:64),[],[],o);
    N = [2:65,70:5:200]';
    x21 = lsqcurvefit(fun1,[0,0,0,0,0,0],N,Hnumerical,[],[],o);
    x22 = lsqcurvefit(fun2,[0,0,0,0,0,0],N,Hnumerical,[],[],o);
    s1 = sum((x11-x21).^2);
    s2 = sum((x12-x22).^2);
    if s1 < s2
        c2 = (c1+c2)/2;
    else
        c1 = (c1+c2)/2;
    end
end